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Abstract 

.-<•  Ur" 

-—''’This  is  the  first  in  a  series  of  three  papers  in  which  discuss  a  method 

<  ^  ' 

for  "post-processing^  a  finite  element  solution  to  obtain  high  accuracy  approxi¬ 
mations  for  displacements,  stresses,  stress  intensity  factors  etc.  Rather  than 
take  the  values  of  these  quantities  "directly"’  from  the  finite  element  solution, 
irk  evaluate  certain  weighted  averages  of  the  solution  over  the  entire  region. 
These  yield  approximations  that  are  of  the  same  order  of  accuracy  as  the  strain 
energy.  AW  obtain  error  estimates,  and  also  present  some  numerical  examples 
to  illustrate  the  practical  effectivity  of  the  technique.  In  the  third  paper 
of  this  series yeT  address  the  matters  of  adaptive  mesh  selection  and 

a 'posteriori  error  estimation. 

/ 


§1  Introduction 

§1.1  The  role  of  post-processing. 

In  many  instances  the  primary  aim  of  a  finite  element  analysis  is  to 
obtain  the  values  of  a  few  important  quantities  with  a  rather  high  accuracy. 
For  instance,  in  structural  mechanics,  the  values  of  displacements,  stresses 
or  stress  intensity  factors  at  a  small  number  of  critical  sites  in  the 
structure  are  important  design  criteria.  Decisions  on  whether  the  structure 
meets  design  specifications,  or  whether  it  is  safe  are  made  on  the  basis  of 
these  few  quantities.  The  bulk  of  the  remaining  numerical  output  of  the 
finite  element  analysis  is  generally  not  scrutinized  so  closely,  but  rather 
is  looked  at  from  a  more  qualitative  viewpoint.  It  may,  for  instance,  be 
used  to  identify  the  critical  points  in  the  construction,  to  obtain  a 
graphical  display  of  the  structure's  deformation  or  to  examine  the  solution's 
plausibility  with  a  view  to  replacing,  if  need  be,  the  particular  mathematical 
model  or  constitutive  law  employed. 

These  considerations  suggest  that  some  thought  should  be  paid  to  how  the 
solution  should  be  "post-processed"  to  obtain  values  for  these  quantities. 
Since  onlv  a  few  quantities  ever  need  to  be  calculated,  we  should  be  willing, 
if  necessary,  to  expend  a  modest  amount  of  computational  effort  on  any  post¬ 
processing  calculation.  A  straightforward  approach  to  post-processing  is  to 
take  the  displacements  or  stresses  directly  as  they  are  output  from  the  finite 
element  computations.  "Curve  fitting"  methods  for  stress  intensity  factor 
calculations  also  fall  into  this  "direct"  category.  There  are,  however,  more 
sophisticated  techniques  which  are  currently  implemented  in  many  commercial 
codes.  We  mention,  for  instance,  the  calculation  of  stresses  at  the  Gaussian 
points  of  certain  element  types,  and  the  use  of  the  "stiffness  derivative" 
method  of  Parks  [  1  ]  or  the  "J-integral"  method  of  Rice  [ 2  1  for  the  deter- 


mination  of  stress  intensity  factors. 

We  shall  show  in  this  and  two  succeeding  papers  that  there  are  post¬ 
processing  procedures  available  which  give  an  accuracy  for  many  physically 
important  quantities  of  the  same  order  as  the  error  in  the  energy  of  the 
finite  element  approximation.  For  example,  in  solving  the  plate  problem 
(that  is,  the  two  dimensional  biharmonic  equation)  with  conforming  elements 
of  degree  p^5  it  is  possible  (with  an  appropriate  finite  element  mesh) 
to  determine  at  any  point  the  displacement,  rotation,  moment  and  shear  force, 
all  with  an  order  of  accuracy  0(N  ^P  1^)  ,  where  N  is  the  number  of  degrees- 
of-freedom  of  the  finite  element  model.  Compare  this  with  the  "direct" 
approach  which  gives  displacements  to  0(11  ^p+^^2),  rotation  to  Q(N  ~^p^) 
and  moments  to  0(N-^P~*^2)  .  The  fact  that  the  "direct"  approach  results  in 
higher  accuracy  for  displacements  than  moments  is  widely  known.  We  remark 
also  that  the  post  processing  approaches  that  we  shall  outline  can  be  well 
understood  from  within  the  standard  "energy"  theory  of  finite  elements.  This 
contrasts  with  the  complex  mathematical  theory  which  is  needed  to  analyze  the 
direct  approach  to  the  determination  of  pointwise  quantities. 

SI. 2  A  general  form  for  post-processing  calculations 

Let  us  denote  by  <t  the  quantity  (e.g.  stress  at  a  point)  which  we  wish 
to  determine.  We  shall  write  $(w)  for  4>  to  indicate  that  4>  relates  to 
a  problem  whose  exact  solution  is  w  .  Now,  let  us  suppose  that  we  know  a 
function  £  so  that 

(1.1)  *  ■  4>(w)  ■  |  w£  dfi  +  R 

n 

where  Q  is  the  region  on  which  our  problem  is  posed,  and  R  is  an  Integral 
which  may  be  computed  using  only  the  problem’s  input  data  (applied  tractions 
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and  forces  etc.).  As  a  simple  illustration  of  what  (1.1)  may  describe, 
consider  the  case  where  $  is  the  value  of  w  at  some  given  point.  Were 
the  influence  function  (Green's  function)  known  for  this  point,  then  we  could 
take  £=0  and  <J>  would  be  expressible  in  terms  of  the  input  data  alone.  So 

(1.1)  would  take  the  form  4'=R  .  Of  course,  the  influence  function  is  not 
in  general  available.  At  the  other  ext reme ,  if  we  take  £  to  be  the  Dirac 
delta  function  at  the  point  under  consideration  then  we  could  set  R=0  and 
<t  =  j  w£  dfl  .  This  is  simply  a  formal  way  of  saying,  evaluate  w  directly 
at  tne  point.  However,  as  we  shall  see  later,  there  are  many  choices  for  £ 
between  these  two  extremes.  Let  us  also  mention  that  (1.1)  may  also  be 
written  in  other  equivalent  forms  (e.g.  after  an  integration  by  parts  we  are 
able  to  obtain  an  integral  in  terms  of  the  derivatives  of  w  instead  of  w 
itself).  The  well  known  J-integral  of  Rice  [2  ]  could  be  thought  of  as  arising 
from  such  a  modified  version  of  (1.1).  The  path  independence  property  of  the 
J-integral  then  corresponds  to  different  choices  for  £  . 

% 

Having  (1.1)  suggests  an  obvious  method  of  approximation  for  $  .  If  w 
is  a  finite  element  approximation  to  w  ,  then  we  could  try  to  approximate 
$  =  <t(w)  by 

(1.2)  $  »  ^(w)  ■  |  w£  dfi  +  R  . 

n 

The  difference  between  4>  and  ^  is  given  by 

(1.3)  e-4>-^  =  J(w-w)£  dft  , 

n 

and  we  see  clearly  that  the  choice  of  £  affects  the  magnitude  of  this  difference. 

%  'C 

If  £  is  the  Dirac  delta  function,  then  $  is  the  point  value  of  w  and  e 
is  the  pointwise  difference  between  w  and  w 
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If  however  £  is  not  concentrated  at  one  point,  then  e  becomes  some  weighted 

average  of  w-w  .  It  is  well  known  that  the  finite  element  method  appears 

more  reliable  when  its  accuracy  is  measured  in  an  average  rather  than  a  point- 

wise  sense.  Spurious  oscillation  which  may  cause  serious  loss  of  accuracy  in 

pointwise  values  of  the  approximate  solution  (especially  of  its  higher 

derivatives)  are  filtered  out  by  averaging.  This  is  especially  the  case  with 

the  p-version  of  the  finite  element  method.  So,  even  at  this  early  stage,  we 

see  that  choices  of  ^'s  which  have  large  support  are  likely  to  give  superior 
<v 

approximations  $  .  A  good  choice  for  £  is  an  important  feature  of  any 
successful  implementation  of  (1.2).  Numerical  experience  however  has  shown 
that,  provided  £  meets  a  few  simple  criteria,  4>  is  quite  insensitive  to  the 
choice  of  £  . 

Another  implementation  issue  that  we  shall  address  is  the  optimal  choice 

% 

of  mesh  (and  consequently  of  w)  for  use  in  (1.2)  This  is  especially  significant 
for  adaptive  finite  element  codes,  where  some  adaptive  criteria  must  be  set. 
Obviously,  if  a  good  approximation  to  $  is  our  ultimate  goal,  this  adaptive 
criteria  should  be  directed  towards  producing  a  w  that  perforins  well  in  (1.2). 
From  a  computational  point  of  view,  evaluation  of  $  needs  at  most  0(N) 

operations,  while  the  solution  of  the  finite  element  problem  itself  usually  needs 

2  7/3 

about  0(N  )  operations  in  two  dimensions  and  Q(N  )  in  three  dimensions. 

Since  only  a  few  evaluations  are  ever  needed,  the  computational  effort  entailed 

is  relatively  insignificant. 

§1.3  Outline  of  the  paper 

This  is  the  first  of  a  series  of  three  papers  which  shall  deal  with  post¬ 
processing  in  the  finite  element  method.  In  this  paper  we  detail  some 
particular  applications  of  the  general  theory  outlined  in  §1.2.  In  §2 


wc  treat  the  example  of  a  one-dimensional,  elastic  supported  beam.  Here  we 
shall  be  interested  in  post-processed  values  for  the  displacement,  rotation, 
moment  and  shear  force  at  certain  points,  along  with  the  average  displacement 
over  a  subsection  of  the  beam.  §3  deals  with  the  simple  two  dimensional  problem 
of  a  membrane  on  an  elastic  support.  For  this  problem  we  shall  be  interested  in 
the  d isplacements  and  stresses  at  certain  points.  Finally  in  §4  we  discuss  a 
numerical  example  related  to  the  problem  treated  in  §3. 

The  second  paper  of  the  series  will  be  devoted  to  applications  in  linear 
facture  mechanics.  In  particular,  we  shall  be  concerned  with  the  computation 
of  stress  intensity  factors  (including  both  and  for  mixed  mode 

fracture).  In  the  third  paper  we  address  the  issues  of  a'posteriori  error 

^  'X, 

estimates  for  4  and  adaptive  mesh  selection  for  w  . 
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§2.__  A  Ono-d imons ional  Kxampjl_e 
§2.1  Formulation  of  the  example 

By  way  of  a  one-dimensional  application  of  the  techniques  outlined  in 
§1.2  we  shall  consider  the  problem  of  a  clamped  beam  (of  unit  length)  on  an 
elastic  support.  The  governing  differential  equation  is 

(2.1)  (aw")"  -  (bw')’  +  cv  =  f  on  (0,1) 

with  the  boundary  conditions 

(2.2)  w(0)  =  w' (0)  =  0  and  w(l)  =  w'(l)  =  0  . 

We  shall  assume  that  a,  b,  c  and  f  are  smooth  functions  which  satisfy 

(2.3)  a(t)  >  a  >  0  (0<t<l) 

b(t) ,  c ( t)  >  0  . 

The  coefficient  a  is  the  rigidity  of  the  beam  while  b  and  c  relate  to 
the  elastic  properties  of  the  support.  For  this  problem  we  shall  be  interested 
in  the  evaluation  of  the  following  five  important  mechanical  quantities: 


(I) 

The 

displacement  of  the  beam  (w)  =  w 

at  0  <  t  <  1  . 

(II) 

The 

rotation  of  the  beam  $2  =  $2(w)  =  w' 

at  0  <  t  <  1 

(III) 

The 

bending  moment  =  $-j(w)  =  aw"  at 

0  <  t  <  1  . 

(IV) 

The 

shear  force  $>£  =  (w)  =  (aw")'  -  bw' 

at  0  <  t  <  1 

(V) 

The 

average  displacement  of  the  subsection 

$5  =  $5(w)  =  (t2-t1)  1  j  2  w  dt 

tl 

jc  t  <_  t2  of  the  beam, 
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§2.2  Expressions  for  ( i  ^  1 ,  ..,4) 

For  the  moment,  let  4>  be  any  function  defined  on  (0,1)  which  satisfies 
the  boundary  conditions  (2.2).  Suppose  also  that  *  is  sufficiently  smooth 
to  allow  any  operations  that  we  carry  out.  Now,  multiply  (2.1)  by  i p  and 
integrate  by  parts  four  times  over  the  entire  interval. 


|  f  <J>  dt  =  j  ((aw")"*  -  (bw')'*  +  cw*)  dt 


[(aw")'*  -  aw"  &  +  aw'*"  -  wCa*")'] 


t-0 


t-0 


_  t+0 

t  1 

[bw'*  -  bw*']  +  (  I  +  |  )  L[*]w  dt 

7+0  l  L 

t 


=  [w(-(a*")'  +  b*')  +  w'(a<f>")  +  aw"(-*')  + 


t-0 


(2.4)  ((aw")'  -  bw')*]_  +  ([  +  f 

t+0  l  jz 


)  L[*]w  dt 


where 


LI*]  =  (a*")"  -  (b*')'  +  c*. 


Let  us  now  be  more  specific  about  the  behavior  of  *  near  t  .  Depending 
upon  which  derivatives  of  *  are  continuous  at  t  ,  (2.4)  will  form  the  basis 
of  our  expressions  for  (i=l,  ..»4)  . 


Case  (I):  Suppose  that 


[  *(1)  ]  I“°  -  0  (i-0,1,2) 


(2.5a) 


t+0 


I  *(3)  1  t_° 

t-o 


-  a(t)"1  , 


while 
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then  (2.4)  gives 


>^(w)  =  w(t)  =  -(  |  +|  )L[$]w  dt  +  |  f<P  dt  . 


0  t 


This  is  exactly  in  the  form  (1.1)  with  R  =  f$  dt  and  t,  =  — L[  4>  ) 


-f 


0 

Notice  that  I.[  d>]  will  in  general  be  discont inuous  at  t  .  We  shall  see 
later  that  from  a  numerical  point  of  view  it  is  important  that  £  be  smooth 
on  (0,1).  So  let  us  append  to  (2.5a)  the  condition 


, t-0 

(2.5b)  [LU]U>]_  =0  j=0,...,n 

t+0 

w’hore  n  is  some  integer,  which,  for  the  moment,  will  remain  arbitrary. 

If  we  select  for  <p  the  influence  function  (Green's  function),  then 
(2.5a)  and  (2.5b)  are  satisfied.  Indeed,  L I <f> ]  =0  on  (0,t)  and  (t,l), 

fl 

and  we  have  <J>  (w)  =  f  $  dt.  Of  course,  in  j^neral,  we  cannot  find  the 
1  J0 

influence  function.  Nevertheless,  functions  <f>  which  satisfy  (2.5a)  and 
(2.5b)  are  readily  constructed,  as  the  following  example  shows: 


Example  1:  We  shall  construct  a  function  $  which  fulfills  all  the  necessary 
conditions,  with  n=l  in  (2.5b).  The  construction  is  done  in  a  number  of 
steps.  First,  define 


Co  o  <  t  <_  t 

\a(t)-1  t  <  t  <  1 


and  then  set 


t  y  x 


♦x(t)  = 


<}>Q(s)ds  dx  dy 


0  0  0 

So  (t)  satisfies  (2.5a).  To  meet  the  requirement  (2.5b)  with  n=l  we 
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Jo  f  i  tie 


,  (t) 


+• 


(t-t)' 


(t-t)' 


0  *  t  <  t 
t  •  t  <  1 


whore  the  coefficients  u  and  6  are  to  he  chosen  so  that  (2.5b)  holds.  It 
is  easy  to  see  that  this  needs 


LlV 


L+0 


-l 


6  =  5  !  a  (T)  (U:l' 


O) 


r«  +  3  a’(t)4!° 


and 


(This  step  can  be  extended  in  an  obvious  fashion  to  handle  n  >  1.) 

We  now  have  a  function  that  satisfies  (2.5),  however  it  may  not  yet 
satisfy  the  boundary  conditions  (2.2).  There  are  many  ways  this  can  be 
remedied.  Two  possibilities  are: 

(i)  Let  x  be  a  smooth  function  which  vanishes,  along  with  its  first 
derivative,  at  t=0  and  t=l  ,  but  has  a  value  of  1  in  an  interval  about 
t  .  Then  <{>  =  meets  all  the  requirements.  We  shall  refer  to  x  as  a 

cut  off  function. 
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(ii)  Let  <  ^  be  a  smooth  function  such  that  <J>j*^(0)  =  and 

$2  ^(l)(i=0,l)  .  Then  we  may  set  $  =  *  Wt'  sball  refer  to  <p^  as  a 

blending  function.  Various  mixtures  of  these  two  techniques  are  possible; 
tor  example,  using  a  cut  .off  function  to  impose  the  boundary  condition  at 
one  end,  and  a  blending  function  approach  to  handle  the  other  endpoint. Q 

The  above  example  typifies  a  general  method  of  construction  that  we  shall 
employ,  either  explicitly  or  implicitly,  throughout  this  series  of  papers: 
Firstly,  a  relatively  simple  function  is  constructed  that  behaves  in  some 
prescribed  "singular"  manner  near  a  given  point.  This  function  must  then  be 
modified  to  ensure  that  it  satisfies  a  set  of  boundary  conditions  on  the 
entire  boundary  of  the  region  of  interest.  Either  cut  off  function,  blending 
function  or  a  combination  of  these  techniques  may  be  used  to  achieve  this.  Let 
us  note  at  this  point,  that  blending  function  techniques  will  usually  lead 
to  a  smoother  modified  function  than  will  the  use  of  a  cutoff  function.  For 
this  reason,  in  a  numerical  setting,  the  blending  function  approach  is  to  be 
preferred. 

Case  II:  Suppose  that 

(i) 

[  *U'  ]_  =0  (i=0,l)  , 

t+0 

1-0  ~  -1 

(2.6a)  [  <j>"  ]_  =  a(t)  1  ,  and 

t-0 

t-0 

i  (ary  ]_  =0 

t+0 


then  (2.4)  gives 
$2(w)  *  v’(t) 


1 

-([  +  j  )  U<f]  w  dt  +  j 
0  t  0 


f  $  dt  . 
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As  for  Case  I,  it  will  turn  out  to  be  important  to  add  the  condition 


(2.6b) 


l  L($](J) 


t-0 

L  -  o 

t+0 


j=0, . 


•  »n  . 


Example  2:  We  shall  construct  a  function  $  that  meets  our  requirements.  We 
proceed  much  as  in  Example  1.  First  we  define  as  in  that  example,  but 

now  set 


4> 


1 


t  X 


0  0 


4>q  (s )  ds  dx  . 


To  satisfy  (2.6b)  (with  n=l,  say)  and  the  third  part  of  (2.6a)  we  may  define 

0  c  t  £  t 
t  <  t  <  1 


*2(0  = 


♦2(t)  +  a(t-t)3  +  6(t-t)4  +  y(t-t)5 


♦x(t) 


where  a,g  and  y  are  to  be  chosen  so  that  (2.6b)  and  the  last  part  of  (2.6a) 
hold.  (That  this  may  always  be  done  follows  since  a(t)^0.)  We  then  proceed 
to  impose  the  boundary  conditions  (2.2)  by  the  same  methods  described  in 
Example  1.1  I 


Cases  III  and  IV:  If  t  is  internal  to  the  interval  then  selecting 

t-0 

f  ♦  ]__  “  0  . 

t+0 

t-0 

(  ♦'  ]_  =  -1  , 
t+0 

t-0 

(2.7)  f  ♦"  ]_  =  0  , 

t+0 

t-0 

'  +  b$']_  «  0  ,  and 

t+0 

[  3_  “0  j-0, . . .  ,n  , 

t+0 
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gJ  VCS 


$^(w)  =  aw"( t )  =  -  ( 


t  1 

M_ 

o  t 


1 


)  L($]w  dt  +  f<p  dt 


with  smooth  L[$l  .  While  choosing  <J>  to  satisfy 


(2.8) 


t-0 

I  <t>  ]_  -  1 

t+0 

7-0 

t  ♦  "  ]_  =  0 

t+0 

7-o 

[  a #" ]_  =  0 

t+0 


t-0 

[  -(a*")'  +  b*' 

]  =  0  ,  and 

t+0 

,2J 

/-N 

-©■ 

c 

* 

« 

o 

II 

•«"*> 

o 

II 

leads  to  the  expression 


4>4(w)  =  ((aw")'  -  bw')(t) 


t  1 

-  ( j^+ j_)L[<f>Jw  dt  + 


dt 


for  <t4(w)  . 

To  treat  these  two  cases  when  t  is  one  of  the  endpoints  we  need  a 
slightly  different  argument.  For  definiteness,  suppose  that  t  is  the 
right  hand  endpoint  t=l  .  We  now  let  $  be  a  smooth  function  on  (0,1) 
which  satisfies  the  boundary  condition  $(0)  *  <f>'(0)  ■  0  .  Analogous  to 
(2.4)  we  now  have 

I  f<p  dt  =  j  ((aw")"<p  -  (bw’)'#  +  c  w$)  dt 
0  0 


12 


(2.9) 


((aw")'  41  -  aw"  4>')it  =  1  +  j  L[ 4>]wdt  . 

0 


If  we  further  set 


(2.10) 


<K1)  =  0  and  4>*  (1)  =  -1 


we  obtain  after  rearranging  (2.9) 

d^(w)  =  aw"(l)  =  “|  Ltd]  wdt  + 


fd  dt  , 


while,  using  in  place  of  (2.10) 


(2.11) 


d(l)  =  1  and  d'(l)  =  0 


leads  to  the  expression 


(w)  =  (aw 


”)’(!)  =  -j1  Ltd]  wdt  +  j1 


Ltd]  wdt  +  j  fd  dt 
0  0 


Functions  d  satisfying  (2.7),  (2.8),  (2.10)  or  (2.11)  and  the  appropriate 
boundary  conditions  are  readily  constructed  using  techniques  similar  to  those 
described  in  Examples  1  and  2. 


§2.3.  An  expression  for  d,. 


The  definition  of  d,.  in  §2.1  is  already  in  the  form  (1.1),  for  we  may 


certainly  write 

(2.12) 


dc  =  dc(w)  =1  w  c  dt 
5  5  i0  ° 


J  n 


where 


(t2"tl) 


-1 


fcl  1  c  1  t2 


otherwise 
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However,  from  our  point  of  view,  this  is  an  unsuitable  expression  for  4>  since 


is  not  smooth.  To  overcome  this  failing  we  may  proceed  as  follows:  Define 


t  s. 


S3  fS 2 


l(t^  (  {  a(s  )  {  {  r  (s  )ds  ds 

r\  J  r\  r\  ”  -*■ 


o  o 


0  0 


2  ds3  ds4 


Let  <f>  be  a  smooth  blending  function  that  satisfies  the  boundary  conditions 


4^(0)  -  ^Ji)(0) 


and 


*21>(1)  =  ^^(l)  (i-0,1) 


and  set  <J>  =  <|>^  -  $  .  Multiply  (2.1)  by  4>  and  integrate  by  parts  four 


times  over  the  entire  interval, 

')"<(>  -  (bw')'<j>  +  c<t>w  dt 


|  f ♦  -  |  (aw") 


0  0 
fl 


L[$]w  dt 


L[<t>2]w)dt 


=  ( (L[<!>.  ]w 

J0 

=  f  U  w  -  (b<t>^')w  +  c*  v  -  L[<jj2)w)dt 

*  A  ^ 


Upon  rearrangement  then 


4>5  =  $5(w)  =  |  ( (b4>^ ' ) '  -  c +  L[(Ji2])wdt  +  j  f  4>  dt 


which  is  of  the  form  discussed  in  §1.2  with 


C  *  (b<f> ,  * )  *  -  cifi,  +  L[$9]  and  R  *  |  f  <p  dt 

0 


r 


In  contrast  to  c  ,  notice  that  £  has  a  continuous  first  derivative  (though  , 


in  general,  a  discontinuous  second  derivative  at  t^  and  t2  ) .  A  t  with 


more  smoothness  can  be  constructed  by  iterating  the  above  process. 


52.4  The  accuracy  of  the  approximations  4\  ( i  =  l , • • » 5) 

In  §2.2  and  §2.3  we  derived  some  integral  expressions  for  the  . 

These  fitted  into  the  general  pattern  discussed  in  §1.2.  We  shall  now  address 
the  important  question  of  the  accuracy  of  the  approximations  =1>^(  w)  which 
arise  when  the  finite  element  solution  w  is  used  in  place  of  w  in  these 
expressions . 

To  be  definite,  suppose  we  have  set  up  a  finite  element  model  of  (2.1)/(2.2) 

using  polynomial  elements  of  degree  p(^3).  Write  I^,...,I  for  the 

intervals  which  comprise  the  finite  element  mesh.  Denote  by  S  the  set  of 

all  admissible  finite  element  functions.  Let  h  =  max  (length  I.  )  .  In  the 

k  k 

case  of  the  problem  (2.1)/ (2. 2)  the  fundamental  orthogonality  property  of  the 

■v 

finite  element  error  w-w  takes  the  specific  form 

(2.13)  [  (a(w-w)"v"  +  b(w-w)'v'  +  c(w-w)v)  dt  =  0 
J0 

for  all  v  in  S  .  Denoting  by  F.(*)  the  energy  expression 

E  (•)  =  I  (a( ( • ) ") 2  +  b ((•)*) 2  +  c  (*)2)  dt 
J0 

the  standard  finite  element  error  estimate  may  be  written  as 

(2.14)  E(w-w)  <_  min  E(w-v*) 

v*€S 

where  the  minimum  is  taken  over  all  v*  from  S  . 

For  the  purposes  of  the  analysis,  we  need  to  introduce  the  auxilliary 
function  which  satisfies 

(2. 15. a)  L[i{f]  =4  on  (0,1)  and 

«K0)  =  *'(0)  -  0  -  <K1)  -  <|',(1>  , 

or  what,  after  integration  by  parts,  is  the  same  thing 
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(2.15b) 


j(a*"  u' 


"  +  b  $i' u'  +  c  u) dt  =  J  C  u  dt 

0 


1 

=  C  u 

4 


for  all  u  with  boundary  values  as  in  (2.2).  Now,  recalling  (1.3),  we  have  the 

% 

following  estimate  for  the  error  $  -  $ 


r  1  ^ 

=  $  -  %  =  C,  (w-w)  dt 


=  |  (a(w-w)"i p"  +  b(w-w)  'tp1  +  c(w-w)(p)  dt 


0 

(putting  u  =  w-w  in  (2.15b) 


■£ 


(a(w-w)"  (ip-v)"  +  b(w-w)'  (ty-v)  '  +  c(w-w)  (4>-v))  dt 


for  any  v  from  S  (by(2.13)).  So 


4>-%l  £  min  ({  (a(w-w)"  (ip-v)"  +  b(w-w)'  +  c(w-w)  (^-v))  dt) 

v€3 


<  min(E  (w-w)^  E  (ij/-v)^) 
~  v€S 


(2.16)  <_  min(E  (w-v*)  2)rain(E(v~v)  5) 

v*€3  s 


using  (2.14).  In  words  then,  the  error  in  <*>  is  bounded  by  the  product  of  the 
energy  norm  difference  between  w  and  its  best  approximation  from  S  , 

and  the  energy  norm  difference  between  and  its  best  approximation  from  S  . 
The  importance  of  the  smoothness  of  C  can  now  be  appreciated.  Smooth  functions 
C  will  give  smooth  auxiliary  functions  ^  ,  and  these  will  be  approximated 

well  by  the  functions  in  S  . 

Let  us  now  try  to  obtain  an  asymtotic  rate  of  convergence  for  ^  which 
will  be  applicable  to  both  the  h  and  p-versions  of  the  finite  element 
method.  First,  recall  the  approximation  result  (see  [3]):  If  z  is  a 
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function  defined  on  (0,1)  and  if  the  smoothness  measuring  quantity 


s  N 
(  l  l 
i=0  k=l 


dt)li 


is  finite  for  some  integer  s  >_  2  ,  then 


(2.17) 


min  E(z-v)  C1  h^”*  ||zj|^ 

v€S  p  2*(s-2)  s 


where  is  a  constant  which  does  not  depend  on  the  function  z  ,  the 

finite  element  mesh  or  the  order  of  elements  used;  and  m=min(p-l,s-2) .  If 

the  load  f  in  (2.1)  and  the  function  C,  are  smooth  enough  to  guarantee  that 

||f||  <  °°  and  |  |s|  |  <  °°  ,  then,  provided  the  coefficients  a,  b  and  c 

S1  S2 

are  sufficiently  smooth,  it  may  be  shown  that 


(2.18) 


M!S]l+4  i  C2  ||f||Si  and 


'*"s2+4i  C3  MCll32 


where  and  do  not  depend  on  f  and  4  respectively.  Having  these 

smoothness  properties  of  w  and  \p  ,  we  may  make  use  of  the  approximation 
bounds  (2.17)  in  (2.16)  to  obtain 


h(W 


(2.19) 


*-*l  i  IMI.^  h*ll,i+t 


(ra.-t-m,) 

icie2c.”«r4?,|i,ii.1 1  mi, 

p  i  z  1 


where  min  (p-l.Sj+2)  (i=l,2)  .  Likewise  (2.14)  gives 

e<„-5)  <  c  c2  ||.||*  . 

pi  1 
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We  mentioned  in  §2.2  and  2.3  that  the  relevant  £'s  could  be  made 
arbitrarily  smooth  (i.e.  n  in  (2.5b)  etc.  was  arbitrary).  For  large 
values  of  n  this  could  become  laborious.  Note  that  (2.19)  allows  dis¬ 
continuities  at  mesh points  without  any  adverse  effects  on  the  accuracy  of 
%  .  What  is  important  is  that  C  be  smooth  in  the  interior  of  each  of  the 
I  .  In  the  h-version  of  the  finite  element  method,  there  would,  at  least 
f rom  our  analys is ,  seem  to  be  no  reason  to  proceed  any  further  than  the  stage 

at  which  s^  +  2  =  p-1  .  At  this  stage  ]<J>-^j  =  0(hml  +  p_1)  ,  and  this  rate 

2m, 

f\j  J_ 

would  not  be  improved  by  increasing  s^  .  Comparing  this  with  E (w-w)  =  0(h  ), 

we  see  ,  as  was  previewed  in  51,  that  the  error  is  at  least  of  the  same 

order  as  the  energy  of  the  error  in  the  finite  element  solution.  For  the  p- 


version  there  would  seem  to  be  no  such  limit.  We  may  increase  s_  indefinitely 


always 


improving  the  convergence  rate  for  4>  as  we  go. 
-(s+s  .,+A)  ^  -2(s  +2) 

0(p  )  and  E  (w-w)  =  0(p  )  . 


Here  we  have 


Of  course,  actual  computations  must  be  carried  out  working  from  only 

a  limited  range  of  non  zero  h's  and  finite  p's  .  In  such  a  setting, 

the  asymtotic  rate  of  convergence  alone  is  not  necessarily  a  good 

indication  of  an  approximation's  accuracy.  As  (2.19)  shows,  the  error 

a, 

in  4  is  related  not  only  to  p  and  h  ,  but  also  to  |U|| 

s2 

and  the  constants  .  As  S£  increases,  the  numerical  values  of  these 

quantities  may  also  increase  dramatically  with  the  net  effect  that  there  is  a 

loss  of  accuracy  in  4>  .  In  practice,  it  is  usually  more  important  to  ensure 

that  the  numerical  value  of  | | £ | |  is  reasonable,  than  to  construct  c's 

s2 

with  high  orders  of  continuity.  Recall  also  our  comment  earlier  that  blending 
function  techniques  are  generally  superior  to  cut-off  function  methods  in  this 
respect.  Another  important  practical  consideration  is  the  choice  of  a  finite 


18 


clement  mesh.  As  (2.16)  shows  the  accuracy  of  $  is  related  to  the  approxima- 
bility  of  both  w  and  i|>  .  So,  an  optimal  mesh  for  calculating  $  would  be 
one  which  was,  in  some  way,  simultaneously  good  for  both  the  original  problem 
(2.1)/ (2. 2)  and  the  auxiliary  problem  (2.15a).  We  shall  explore  this  question 
in  the  third  of  this  series  of  papers.  The  sorts  of  concerns  touched  upon  in 
this  paragraph  are  of  great  importance  in  two  dimensional  problems.  For  many 
such  problems  there  is  no  complete  analogue  to  2. 18),  and  we  are  denied  the 
luxury  of  being  able  to  make  \p  as  smooth  as  we  wish. 
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§3  T wo  D linens  i o nal  Problems 
§3.1  Formulation  of  the  problem 

To  illustrate  the  ideas  of  §1.2  in  a  two-dimensional  setting,  we  shall 
consider  in  some  detail  the  simple  model  example  of 

2 

(3.1)  v  w  -  kw  »  f  in  SI  ,  a  polygonal  region, 

with  the  boundary  condition 


(3.2)  w-0  on  3  0  ,  the  boundary  of  Q  . 

Here  we  suppose  that  k  ^  0  and  assume,  for  simplicity,  that  k  is  a 
constant  and  f  a  smooth  function.  The  problem  (3.1)/(3.2)  could,  for 
example,  be  thought  of  as  describing  a  polygonal  membrane  on  an  elastic  support, 
that  is  fixed  along  its  edges.  We  shall  be  concerned  with  evaluating  the 
following  quantities  which  are  related  to  w  : 

(i)  The  displacement  4^  =  4>^(w)  =  w(x)  at  a  point  x  =  (x^,x2)  in  ft  . 

(ii)  The  stress  $2  =  <?2  (w)  =  Vw-fl(y)  at  a  point  y  =  (y^y^  on  3P  ,  which 
is  "far"  from  a  corner  point.  (The  case  of  y”  "close"  to  a  comer 
point  will  be  discussed  in  our  second  paper.) 


§3.2  An  integral  for  <J>^ 

Let  <)>(x)  be  an  arbitrary  function  defined  and  sufficiently  smooth  on 
P  -  {x}  ,  which  vanishes  on  3P  .  For  c  >  0,  small  enough,  denote  by  Sc 
a  disc  with  centre  x  and  radius  e  which  lies  in  P  .  Multiply  (3.1)  by 
<p  and  integrate  over  P-S£.  Using  Green's  Theorem,  we  obtain 


(3.3) 


dA 


where  L  [•  3  ■  V^(  •) 

pointing  towards  x 


(Vwn<j>  -  V(J>«n  w)  ds  + 


3S 


L[4]w  dA 


P-S 


c 

k(*)  and  fi  denotes  the  unit  normal  on 
Now,  impose  the  extra  conditions 


3S 


e 
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(3.4) 


where  r(x) 
yields 


<Kx)  =  (2ti)  1log  r(x)  +  0(1) 

Vifi(x)  =  V((2tt)  'log  r(x))+  o(r(x)  ')  j 


-  2  —  2  4 

((x^-Xi)  +  (x2-x  )  )  •  Then,  in  the  limit 


(as  x  -*  x) 


as  c  ->  0  ,  (3. 3) 


(3.5)  4j(w)  =  w(x)  = 


1 .  ( f  ]  w  dA  + 

4 

(1 


ft 

n 


dA 


Note  that  the  integrals  appearing  on  the  right  hand  side  of  (3.5)  are  possibly 
improper.  U’e  see  that  (3.5)  is  precisely  in  the  form  required  by  (1.1)  with 
i ;  =  — L 1 4-  ]  and  R  =  j"  f  <£  dA  .  Notice  also  that  were  <}>  the  influence 

function  (Green's  function)  for  (3.1)/(3.2),  then  the  first  integral  on  the 
right  hand  side  of  (3.5)  vanishes.  In  general,  of  course,  the  influence 
function  is  not  available. 

Just  as  in  §2,  it  will  turn  out  that  from  a  numerical  viewpoint,  it  is 
important  for  t,  to  be  a  smooth  function.  The  problem  of  selecting  a  suitable 
5  ,  or  what  is  the  same  thing,  of  choosing  <f  appropriately  can  be  thought  of 
as  having  two  aspects.  Firstly,  ensuring  that  L[ip)  is  smooth  in  the  immediate 
neighborhood  of  x  ;  and  secondly,  of  imposing  the  boundary  conditions  on  <f> 
in  such  a  way  that  no  unsmooth  behavior  of  t>  is  introduced. 

Let  us  talk  in  more  detail  about  these  points  as  they  relate  to  our  model 
problem.  It  Is  easy  to  verify  that  if 

(3.6)  <£(x)  =  (2u)~1(l  +  ~  r(x)2)  log  r(x)  , 

-v  — 2  —  —  'v 

then  I.[$]  =0(r  log  r)  in  the  vicinity  of  x.  Now  $  has  the 

required  asymtotic  behaviour  (3.4),  however,  it  does  not  vanish  on  8ft  ,  and 

so  cr.nnot  be  used  directly  for  4  in  (3.5).  Let  us  suppose  for  a  moment  that 


x 


is  not  too  close  to  Oft  .  Then,  as  for  the  one-dimensional  case,  there  are 


a  number  of  techniques  for  modifying 


*  .  We  could,  for  instance,  proceed  in 

one  of  the  following  ways: 

(i)  Let  X(x)  he  a  smooth  "cut  off"  function  which  vanishes  on  3f l  , 

rV 

but  is  a  constant,  equal  to  1,  in  a  neighborhood  of  x  .  Then  y(x)  =  X(x)  1  ( >:) 
satisfies  all  our  requirements. 

(ii)  Let  $  (x)  be  a  smooth  "blending"  function  on  :!  which  agrees  with 

'V, 

C  ( x )  on  j£  .  We  could  then  take  4(x)  =  4(x)-<,*  (x)  . 

o 

(iii)  L’s e  a  combination  of  the  above  two  techniques — "cut  off"  functions 
to  handle  part  of  the  boundary,  "blending"  functions  for  the  remainder. 

In  the  case  that  x  is  very  close  to  3.°.  ,  technique  (i)  is  not  the  proper  method 
to  use,  an  the  "cut  off"  function  X  would  then  have  large  derivatives 

in  the  region  between  x  and  3S1  .  Methods  (ii)/ (iii)  provide  better  ways  of 

■v 

handling  this  case.  Note  however,  that  though  I.[f]  is  smooth  near  x,  4 

is  not.  For  an  arbitrary  "blending"  function  < ,  agreeing  with  4  on  3Q 

near  x  .there  is  no  reason  to  expect  that  LI  4q]  be  smooth.  One  way  around 

this,  at  least  in  the  case  when  x  ,  though  close  to  the  boundary  of  W  ,  is 

% 

far  from  a  corner  point  of  the  boundary,  is  to  formally  extend  4  across  the 
straight  line  segment  of  the  boundary  closest  x  .  Now  let  ip  be  the  reflection 
of  this  extension  back  into  Q  .  Then  L  ip  is  smooth,  and  4  agrees  with 
4  on  the  straight  line  segment  of  the  boundary  closest  x  .  Standard  "blending" 
techniques  may  then  be  used  to  deal  with  the  remaining  three  sides  of  Q  . 

There  are  obvious  extensions  of  the  above  ideas  to  domains  with  curved  boundaries. 

§3.3  An  Integral  for 

For  definiteness,  suppose  that  y  lies  on  the  straight  line  segment  ((1,-1), 

(1,1))  which  forms  part  of  311  .  This  time,  let  if  be  an  arbitrary,  sufficiently 

smooth  function  defined  on  H  ,  which  vanishes  on  3fl-{x)  .  For  e  >  0  »  small 

enough,  denote  by  S+(x)  a  half  disc  with  centre  x  and  radius  e  ,  which  lies 

c 

in  fi  .  Multiply  (3.1)  by  41  and  Integrate  over  fi-S*  •  Using  Green's  Theorem 
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wo  obtain 


(3.7) 


f  <  dA  = 


(Vwiir  * nvs){js  + 


L[$]  w  dA  , 


fi-S 


Q-S 


whore  denotes  the  circular  portion  of  the  boundary  of  ,  and  as  usual 


n  denotes  a  unit  normal  pointing  towards  the  centre  of  S^_  (see  Fig.  1).  No 


i f  we  impose  the  extra  conditions 


(3.8)  ,(x)  =  A -c.°iL— +  o(r(x)“  ) 

*  rfx) 


(as  x  -*•  y 
from  v;ithin  Q  ) 


v\  =  1  cos  P(x)  ~  -2. 

v,x'  v(  - — - )  +  o(r(x)  ) 

71  r(x) 


where  r(x),  0(x)  arc  plane-polar  coordinates  centered  on  y  (see  Fig.  1), 
then,  in  the  limit  as  e  -*■  0  ,  (3.7)  gives 


(3-9) 


V“>  "  £</>  - 


I 


( 


U*J  W  dA  -  f  {  dA  . 


This  is  in  the  form  of  (1.1)  with  C  =  L[$]  and  R  =  -jf  4>  dA  .  In  general 


the  integrals  on  the  right  hand  side  of  (3.9)  will  be  improper. 
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As  usual  we  should  choose  4  such  that  £  is  smooth.  This  problem  can 
he  approached  in  an  analogous  fashion  to  that  outlined  in  §3.2.  In  our  case 
it  can  he  verified  that  if 

(3.10)  *<x)  -  (£  U  •-  |  r(x)2)  log  7(x))  , 

then  L [  i,s]  is  smooth  about  x  (in  fact,  L[  4}  =  0  (r(x)  log  r(x))  . 

Mow  ^  has  the  necessary  asymtotic  and  boundary  behavior  near  y  ,  but  it  does 

% 

not  vanish  everywhere  on  3 ft  .  To  construct  suitable  $  s  based  on  $  we  may 
use  obvious  adaptations  of  the  "cut  off"  or  "blending"  techniques  outlined  in 
§3.2. 
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§3.4  More  general  problems  than  (3.1)/(3.2) 

The  methods  of  §3.2  and  §3.3  may  be  applied  to  problems  more  general 
than  (3.3)/(3.2)  (e.g.  problems  in  elasticity,  problems  with  non-constant 
coefficients  or  non-homogeneous  boundary  conditions).  We  shall  not  go 
into  details  here.  Let  us  however,  just  mention  that  for  the  problem 
(3.1)  with  the  essential  boundary  condition  (3.2)  replaced  by  the  natural 
boundary  condition  Vw  •  n  *  0  ,  there  would  have  been  no  need  in  §3.2  to 
impose  any  boundary  condition  on  $  .  The  expression  (3.5)  in  this  case 
would  have  to  include  an  additional  term,  namely,  a  line  integral  around 

an  . 


§3.5  The  accuracy  of  the  approximations 

Suppose  that  we  set  up  a  finite  element  model  of  the  problem  (3.1)/(3.2). 
In  the  usual  way  let  us  partition  ^  into  elements  E^,  E2»...,En  say  (we 
do  not  need  to  be  specific  about  the  shapes  of  the  elements) ,  and  assume  that 

r\, 

on  each  element  we  represent  w  by  a  polynomial  of  degree  p  .  The  demand 
of  conformity  requires  that  these  polynomials  be  continuous  across  the  inter¬ 
element  boundaries,  and  vanish  on  3 .  Denote  by  S  the  set  of  all  such 
finite  element  functions.  Let  be  a  characteristic  linear  dimension 

of  E.  ,  and  set  h  =  max  h.  . 

j  J  j 

'V 

The  finite  element  solution  w  satisfies 


(3.11)  J  (  Vw  Vv  +  k  w  v)  dA  *  -j  f  v 


dA 


for  all  v  from  S  ;  in  addition,  we  have 
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(3.12)  j (V(w-w)  Vv  +  k(w-w)v)  dA  =  0 

for  all  finite  element  functions  v  in  S  .  Defining  the  strain  energy  expression 
by 


E(-) 


(V(*)2  +  k(-)2)dA  , 

n 


we  have 


(3.13)  E(w-w)  <_  min  E  (w-v) 

v  €  S 


In  line  with  the  general  procedure  outlined  in  §1  we  consider  approxima¬ 
tions  \  =  *  (w)  and  (w)  to  and  4>2  .  In  either  case,  we  make 

an  error  of  the  form 


(3.14)  e  =  $  -  =  j  £  (w-w)  dA 

n 

Now,  just  as  for  the  one-dimensional  case,  introduce  an  auxiliary  function 
iji(x)  which  satisfies 

V2i>  -  kij>  *  -c 
i^=0  on  3n  , 

or  equivalently. 


i  (V  i(i  7u  +  k  u)dA 


|  4  u  dA 
Q 


for  all  functions  u  which  vanish  on  3f2  and  for  which  E(u)  is  finite. 

<v 

We  may  certainly  choose  u  =  w-w  to  obtain  from  (3.14) 


*\j 

iji  V(w-w) 


+  k  t{i(w-w)  dA 


> 


1 


e 


and  using  (3.12)  we  see  that  for  any  finite  element  function  v  from  S 

|t'i  =  |  (V  ((/-v)  V  (w-w)  +  k  (£-v)  (w-w))  dA  | 

9. 

L  ^  1, 

<_  i.(^-v)  '  F.  (w-w)  2 

So  on  clioosing  v  to  minimize  E  ( v)  ,  and  recalling  (3.  13),  we  have 

i,  *  i, 

(3.15)  1 1* '  min  F(y-v)  2  min  e(w-v 

v€3  v*  €3 

This  estimate  is  telling  us,  exactly  as  did  (2.16)  in  the  one-dimensional  case, 

'V 

that  the  accuracy  of  i>  depends  on  how  well  both  the  solution  w  of  the 
original  problem,  and  the  solution  i)>  of  the  auxiliary  problem,  can  be 
approximated  in  the  energy  norm  by  the  finite  element  functions  in  S  . 

If  we  try  to  obtain  rates  of  convergence  for  $  ,  we  come  up  against  some 
important  differences  between  the  one  and  two  dimensional  cases.  In  general, 
the  analog  of  (2.18)  holds  only  if  the  boundary  of  SI  is  smooth.  If  9fl  is 
not  smooth  then  (2.18)  must  be  modified  to  account  for  some  special  singular 
terms  that  arise  because  of  corners  of  9ft  (see  [  4  ]).  These  singular  terms 
govern  the  smoothness  and  approximability  of  w  and  ip  .  The  analog  of  estimate 
(2.17)  is  also  more  complicated  in  the  two  dimensional  case  (see  [  3]).  None¬ 
theless,  if  the  mesh  has  the  proper  level  of  refinement  around  the  corners  of 
9ft  ,  then  similar  results  to  those  in  the  one-dimensional  case  can  be 
achieved  if  the  rate  of  convergence  is  now  measured  with  respect  to  the 
number  of  degrees-of-f reedom  rather  than  p  and  h  .  To  go  into  further 
details  is  beyond  the  scope  of  this  paper. 


(  *****  .  . 


§ 4  A  Numerical  Example 
14.1  Formulation  of  t lie  example 

As  a  practical  demonstration  of  the  methods  discussed  in  §3,  we  shall 
consider  some  numerical  results  for  the  problem  modelling  a  square,  uniformly 
loaded  membrane  which  is  fixed  along  its  edges.  More  specifically,  we  deal 
with  the  problem  governed  by  the  differential  equation 


(4.1) 


V2w  =-l  on  ft  =  < -1 , 1) 2 


and  boundary  conditions 


(4.2) 


w  =  0  on  the  boundary  3ft  of  ft 


We  shall  employ  the  theory  of  §3  for  the  calculation  of  approximate 
values  for: 

(I)  The  displacement  at  the  centre  of  the  membrane:  =  <f>^(w)  =  w(0) 

(II)  The  stress  at  the  point  P(1,0):  <J>2  =  <t>2(w)  =  ^  (P)  • 

By  the  method  of  separation  of  variables,  an  infinite  series  representation  of 
w  can  be  found.  Using  this  series  the  following  exact  values  (accurate  to 
5  significant  figures)  can  be  calculated: 


dA  =  .56231 


$^(w)  =  w(Q)  =  -.29469 


$2 (w)  -f~(P)  -  -67528 


Let  us  also  note  that  the  solution  w  is  relatively  smooth  (in  fact,  it  has 
square  integrable  second  derivatives,  though  not  square  integrable  third 
derivatives) . 

§4.2  The  finite  element  approximation 

We  shall  consider  a  simple  finite  element  model  of  (4.1)/(4,2).  Namely, 
bilinear  elements  on  a  square  uniform  mesh.  By  the  symmetry  of  the  problem, 
we  need  only  actually  calculate  using  the  quarter-segment  OQRP  of  ft  (see 


0  (0,1) 

RCi.l) 

P  (1 ,0) 


(1,-1) 

Figure  2;  The  region  of  the  model  problem. 

Figure  2).  For  this  problem  we  expect  the  following  rates  of  convergence 

2  'v, 

0(h)  for  the  energy  E  (w) 

0(h)  for  the  energy  norm 

2  'x# 

0(h  )  for  the  displacement  v(0) 

3w 

0(h)  for  the  stress  t —  (P) 

9X1 

where,  as  usual,  h  denotes  the  length  of  the  side  of  an  element. 

Using  a  uniform  mesh  and  elements  of  degree  2  or  higher  we  would  obtain 

3 

0(h  )  for  the  energy 

0(h^'^~E)  for  the  displacement  w(0) 

and  0(hA*  e)  for  the  stress  r — (P) 

9xl 

where  e  >  0  is  an  arbitrary  small  number.  For  the  h-p-version,  it  is 
possible  by  suitable  refinement  about  the  corner  points  to  obtain 
arbitrarily  large  orders  of  convergence  with  respect  to  the  number  of 
degrees-of-f reedom  (see  [3]). 
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§4.3  Calculation  of  ^(w) 

In  accord  with  the  theory  developed  in  §3.2,  we  use  the  formula 

%  %  (  ?  'V/  f 

lf(w)  =  -jV<fwdA  -  J  dA  , 

Si  fl 

where  <f  takes  the  generic  form 


41  =  X(x1>x2)  log  (x^  +  x22)5  -  ^o(x1,x2>] 


We  consider  two  choices  for  : 


Case  (a):  Xix^.x,^)  =  X(x^)X(x2)  where 

f1 

x(t;  =  < 

\i  -  8( | t | -  i)3 

(see  Figure  2  ,  and 


o  5  1 1|  _< 

%  <  I t|  £  1 


Figure  3:  Cutoff  function  used  in  the  evaluation  of  4^ 


case  (a) 


Case  (b):  X(x2,x2)  =  1  ,  and 


VV5^ 


=  27  (l°B 


(l+xj  K1+x22  )\H 


)) 


In  Case  (a)  we  have  employed  a  cut  off  function  technique  to  enforce  the 

boundary  conditions  on  <p  »  while  in  Case  (h)  a  blending  function  method  has 

been  used.  The  first  integral  in  the  for  ’going  formula  for  may  be 

calculated  by  numerical  quadrature.  (We  used  Gaussian  quadrature.)  The 

second  integrand  is  singular  at  0  .  Th  os  integral  may  be  evaluated 

analytically.  However,  it  is  also  possible  to  calculate  it  numerically  by 

2  2  2 

the  following  procedure:  Choosing  p  such  that  V  p  =  1  (e.g.  p  =  +x2  )), 

integration  by  parts  gives 


f  2 

f 

► 

4>  V  p  dA  = 

1)1  Vp*A  ds  - 

p  ds  + 

SI 


an 


an 


n 


p  dA  +  p (0,0) 


All  the  integrals  on  the  right  hand  side  are  nonsingular,  and  may  be  readily 
evaluated  by  numerical  means. 

The  results  of  the  computations  are  shown  in  the  middle  section  of  Table  1. 

>v 

For  comparison,  we  also  list  the  value  of  the  finite  element  solution  w  at  0. 
Notice  that  w(0)  and  both  cases  of  all  show  an  0(h  )  rate  of 

convergence.  This  is  as  we  would  expect.  Observe  also  the  superiority  of 
the  post-processed  value  in  Case  (b)  over  Case  (a).  This  is  in  line  with 
our  previous  comment  that  blending  function  techniques  can  usually  be  expected 
to  perform  better  than  cut  off  function  methods.  (Looking  at  the  definition  of 
X  in  Case  (a)  and  examining  Figure  3  si-  -ws  that  indeed  X  changes  quite 
rapidly  in  the  region  ) X2 I  >  *S  •  In  terms  of  the  arguments  we  presented  in 
§3.5,  we  should  therefore  not  expect  the  corresponding  to  bo  as  well 
approximated  from  within  our  finite  element  subspsce  as  it  would  be  in  Case  (b) — 
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No.  of  elements  in  quarter 

segment  (uniform  mesh). 

4 

16 

64 

Exact 

Va  1  ue 

% 

I.nergy  norm  error  in  w 

-  A-:  (>’-«!  V 

30.1% 

15.2% 

7.62% 

\  E(w)  / 

w(0) 

(*). 310714 

. 298393 

.295596 

(relative  %  error) 

(5.4%) 

(1-3%) 

(.31%) 

(a) 

.268783 

.287205 

.292751 

.29469 

tx(w) 

(8.8%) 

(2.5%) 

(.65%) 

(relative  %  error) 

(b) 

.287306 

.292829 

.294220 

(2.5%) 

(.63%) 

( . 16%) 

- . J 

% 

!>> 

.482142 

.565480 

.616687 

(relative  %  error) 

(29%) 

(16%) 

(8.7%) 

(a) 

.64758 

.67197 

.67463 

(4.1%) 

(.49%) 

(.096%) 

.67528 

(b) 

.66623 

.67313 

.67477 

(relative  %  error) 

(1.3%) 

(.32%) 

(.076%) 

(c) 

.66482 

.67276 

.67468 

(1.5%) 

(.37%) 

(.089%) 

...  ... 

. 

i  .   ] 

TABLE  1.  Table  of  the  results  of  the  numerical  calculations. 

((*):  negative  signs  have  been  supressed  in  this  table) 
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see  (3.13)). 


s,  \ 

The  tact  that  t  Ik1  accurac  ies  of  w(0)  and  S^(w)  art*  comparable  in  this 
example  is  a  eonseipioncr*  of  our  u*inp  bilinear  elements.  Nonetheless,  Table  1 
shows  that  in  Case  (b),  the  values  are  twice  as  accurate-  as  the  w(0) 

values  for  the  same  number  of  elements.  Tutting  this  another  way,  for  the 
same*  accuracy  the  "direct"  displacement  method  would  require*  twice  as  many 
el  enter  t  s  as  the  pest -process  iiig  approach  of  Case  (b). 

54.3  Calculation  of  ‘  (w) 

In  the  case  of  our  model  problem  the  theory  of  §3.3  leads  to 


s,  *v  f  o  *v  [ 

$2  (w)  =  j  Vi  w  dA  +  I  6  dA 


where  4  takes  the  Generic  form 


t  =  X(x1,x2)  (  ) 


xr! 


2  2 

^(Xj-1)  +x2 


^o(xl,X2) 


We  shall  treat  three  cases: 


Case  (a) : 


XCx^x.,)  =  -  Sx^3  +  Sx.^2) 


(See  Figure  4(1)) 


*o(xl’X2) 


xrl 


(Xj-1)  +  1 


Case  (b):  X(x 


1,X2>  ’{-I.,! 


(See  Figure  4 ( i  1 ) ) 


,  .--IC1  — 

0  (xrl)2  +  1 


!i  1  Xj  1  1 

o  £  <  h 

-1  <_  Xj_  <  0 


0  £x^  <_  1 
-1  <  xx  <  0 
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r (cl :  X(xj,x  )  -  1 


x  -i 

,  -  +  (x  -1) 

Uj-1)“  +  1  \ 


4+x 


1  _  1 

-  2-  5 


01  -101 

(i)  (ii) 

'X, 

Figure  A:  Cut  off  functions  used  in  the  evaluation  of 

(i)  case  (a) 

(ii)  case  (b) 


Cases  (a)  and  (b)  correspond  to  our  using  a  blending  function  technique  to 
satisfv  the  boundarv  condition  on  the  edges  x2  ~  i-1  »  and  a  cut  off  function 
method  to  handle  the  edge  *1  =  -1  .  In  Case  (c),  a  blending  function  method 
is  used  to  handle  the  entire  boundary.  Concerning  the  actual  evaluation  of 
^2 ( w)  ,  the  same  comments  made  in  54.2  about  i ^ (w)  apply  here  also. 


The  results  of  the  calculations  are  shown  in  the  lower  part  of  Table  1, 

where,  for  comparison,  we  have  also  listed  the  corresponding  values  of 
-v 

(P)  .  In  contrast  to  the  situation  for  the  displacements,  we  see  that 

3*1 

the  post-processed  values  for  $2  are  markedly  more  accurate  than  the  "direct' 
value  ~~(P)  .  We  see,  as  theory  predicts,  an  0(h)  rate  of  convergence 

(JX. 

-*•  2.  r\j 

the  "direct"  value,  but  an  0(h)  rate  for  $2  • 


for 
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